home *** CD-ROM | disk | FTP | other *** search
/ Enter 2004 January / enter-2004-01.iso / files / maxima-5.9.0.exe / {app} / share / maxima / 5.9.0 / src / rombrg.lisp < prev    next >
Encoding:
Text File  |  2003-02-09  |  4.9 KB  |  150 lines

  1. ;;; -*-  Mode: Lisp; Package: Maxima; Syntax: Common-Lisp; Base: 10 -*- ;;;;
  2. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  3. ;;;     The data in this file contains enhancments.                    ;;;;;
  4. ;;;                                                                    ;;;;;
  5. ;;;  Copyright (c) 1984,1987 by William Schelter,University of Texas   ;;;;;
  6. ;;;     All rights reserved                                            ;;;;;
  7. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  8. ;;;     (c) Copyright 1982 Massachusetts Institute of Technology         ;;;
  9. ;;; Original code by CFFK.  Modified to interface correctly with TRANSL  ;;;
  10. ;;; and the rest of macsyma by GJC                                       ;;;
  11. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  12.  
  13. (in-package "MAXIMA")
  14. (macsyma-module rombrg)
  15. (load-macsyma-macros transm numerm)
  16.  
  17. (declare-top(special user-timesofar))
  18.  
  19. ;;; the following code if for historical frame of reference.
  20. ;;;(defun fmeval3 (x1) 
  21. ;;;       (cond ((integerp (setq x1 (meval x1))) (float x1))
  22. ;;;         ((floatp x1) x1)
  23. ;;;         (t (displa x1) (error '|not floating point|))))
  24. ;;;
  25. ;;;(defun qeval3 (y1 x1 z)
  26. ;;;       (cond (x1 (fmeval3 (list '($ev) y1 (list '(mequal) x1 z) '$numer)))
  27. ;;;         (t (funcall y1 z))))
  28.  
  29. (DEFMVAR  $ROMBERGIT 11. "the maximum number of iterations" FIXNUM)
  30. (DEFMVAR  $ROMBERGMIN 0. "the minimum number of iterations" FIXNUM)
  31. (DEFMVAR  $ROMBERGTOL 1.0e-4 "the relative tolerance of error" FLONUM)    
  32. (DEFMVAR  $ROMBERGABS 0.0 "the absolute tolerance of error"  FLONUM)
  33. (DEFMVAR  $ROMBERGIT_USED 0 "the number of iterations actually used." FIXNUM)
  34.  
  35. (DEFVAR ROMB-PRINT NIL ); " For ^]"
  36.  
  37.  
  38. (defun $ROMBERG_SUBR (FUNCTION LEFT RIGHT)
  39.   (BIND-TRAMP1$
  40.    F FUNCTION
  41.    (LET ((A (FLOAT LEFT))
  42.      (B (FLOAT RIGHT))
  43.      (X 0.0)
  44.      (TT (make-array  $rombergit :element-type 'double-float))
  45.      (RR  (make-array  $rombergit :element-type 'double-float))
  46.      (USER-TIMESOFAR (cons 'romb-timesofar user-timesofar))
  47.      (ROMB-PRINT NIL))
  48.      (setq X (-$ B A))
  49.      (SETF (AREF$ TT 0)
  50.        (*$ x (+$ (FCALL$ F b) (FCALL$ F a)) 0.5))
  51.      (SETF    (AREF$ RR 0.)
  52.         (*$ x (FCALL$ F (*$ (+$ b a) 0.5))))
  53.      (do ((l 1. (f1+ l)) (m 4. (f* m 2.)) (y 0.0) (z 0.0) (cerr 0.0))
  54.      ((= l $rombergit)
  55.       (MERROR "ROMBERG failed to converge"))
  56.        (DECLARE (FLONUM Y Z CERR)
  57.         (FIXNUM L M))
  58.        (setq y (float m) z (//$ x y))
  59.        (SETF (AREF$ TT L) (*$ (+$ (AREF$ tt (f1- l))
  60.                   (AREF$ rr (f1- l))) 0.5))
  61.        (SETF (AREF$ RR L) 0.0)
  62.        (do ((i 1. (f+ i 2.)))
  63.        ((> i m))
  64.      (COND (ROMB-PRINT
  65.         (SETQ ROMB-PRINT NIL)            ;^] magic.
  66.         (MTELL "Romberg: ~A iterations; last error =~A;~
  67.                 calculating F(~A)."
  68.                I
  69.                CERR
  70.                (+$ (*$ z (float i)) a))))
  71.      (SETF (AREF$ RR L) (+$ (FCALL$ F (+$ (*$ z (float i)) a))
  72.                 (AREF$ rr l))))
  73.        (SETF (AREF$ RR L) (*$ z (AREF$ rr l) 2.0))
  74.        (setq y 0.0)
  75.        (do ((k l (f1- k))) ((= k 0.))
  76.      (DECLARE (FIXNUM K))
  77.      (setq y (+$ (*$ y 4.0) 3.0))
  78.      (SETF (AREF$  TT (f1- K))
  79.            (+$ (//$ (-$ (AREF$ tt k)
  80.                 (AREF$ tt (f1- k))) y)
  81.            (AREF$ tt k)))
  82.      (SETF (AREF$ RR (f1- K))
  83.            (+$ (//$ (-$ (AREF$ rr k)
  84.                 (AREF$ rr (f1- k))) y)
  85.            (AREF$ rr k))))
  86.        (setq y (*$ (+$ (AREF$ tt 0.)
  87.                (AREF$ rr 0.)) 0.5))
  88.        ;;; this is the WIN condition test.
  89.        (cond ((and
  90.            (or (not
  91.             (< $rombergabs
  92.                (setq cerr
  93.                  (abs (-$ (AREF$ tt 0.)
  94.                       (AREF$ rr 0.))))))
  95.            (not (< $rombergtol
  96.                ;; cerr = "calculated error"; used for ^]
  97.                (setq cerr (//$ cerr
  98.                        (cond ((= y 0.0) 1.0)
  99.                          (t (abs y))))))))
  100.            (> l $rombergmin))
  101.           (SETQ $ROMBERGIT_USED L)
  102.           #+maclisp
  103.           (progn (*rearray tt) (*rearray rr))
  104.           (return y)))))))
  105.  
  106.  
  107.  
  108. (defun romb-timesofar () (setq romb-print t))  ;^] function.
  109. ;;; Making the ^] scheme work through this special variable makes
  110. ;;; it possible to avoid various timing screws and having to have
  111. ;;; special variables for communication between the interrupt and MP
  112. ;;; function.  On the other hand, it may make it more difficult to
  113. ;;; have multiple reports (double integrals etc.).
  114.  
  115.  
  116. ;;; TRANSL SUPPORT.
  117.  
  118. (DEFPROP $ROMBERG_SUBR $FLOAT FUNCTION-MODE)
  119.  
  120. (DEFUN ROMBERG-MACRO (FORM TRANSLATEP)
  121.   (SETQ FORM (CDR FORM))
  122.   (COND ((= (LENGTH FORM) 3)
  123.      (COND (TRANSLATEP
  124.         `(($ROMBERG_SUBR) ,@FORM))
  125.            (T
  126.         `((MPROG) ((MLIST) ((MSETQ) $NUMER T) ((MSETQ) $%ENUMER T))
  127.               (($ROMBERG_SUBR) ,@FORM)))))
  128.     ((= (LENGTH FORM) 4)
  129.      (LET (((EXP VAR . BNDS) FORM))
  130.        (COND (TRANSLATEP
  131.           `(($ROMBERG_SUBR)
  132.             ((LAMBDA-I) ((MLIST) ,VAR)
  133.                 (($MODEDECLARE) ,VAR $FLOAT)
  134.                 ,EXP)
  135.             ,@BNDS))
  136.          (T
  137.           `((MPROG) ((MLIST) ((MSETQ) $NUMER T) ((MSETQ) $%ENUMER T))
  138.                 (($ROMBERG_SUBR)
  139.                  ((LAMBDA) ((MLIST) ,VAR) ,EXP)
  140.                  ,@BNDS))))))
  141.     (T
  142.      (WNA-ERR '$ROMBERG))))
  143.  
  144. (DEFMSPEC $ROMBERG (FORM)
  145.   (MEVAL (ROMBERG-MACRO FORM NIL)))
  146.  
  147. (def-translate-property $ROMBERG (FORM)
  148.   (LET (($TR_NUMER T))
  149.     (TRANSLATE (ROMBERG-MACRO FORM T))))
  150.